##------------------------------------------##
## Lauren Peritz
## ISQ 2019 - WTO Compliance
## Replication Files Parts 1, 2, and 3                
## Update 2019.10.30
##------------------------------------------##

rm(list=ls(all=TRUE))

# Packages ----------------------------------------------------------------
require(readstata13)
require(ggplot2)
require(stargazer)
require(sampleSelection)
require(sandwich)
require(lm.test)
require(MASS)

# #require(AER)
# #require(effects)
# require(clusterSEs)
# #require(gplot)
# require(rms)
# require(effects)
# #require(xlsx)
# require(xtable)


## Data Input --------------------------------------------------------------

setwd("~/Dropbox/WTO_Compliance")

# Main data set for Tables 1 and 2 and robustness checks
d1 <- read.dta13("Replication/Peritz_ISQ2019_WTO1.dta")

# Data for Heckman Selection robustness checks
d2 <-read.dta13("Replication/Peritz_ISQ2019_WTO2.dta") 

# Data from Davis (2012) for Selection into Litigation
d.davis <- read.dta13("Replication/Davis2012_replication_progress.dta")

# Data for single-step linear regression in robustness checks
d3 <- read.dta13("Replication/Peritz_ISQ2019_WTO3.dta")


######################################################################################
#### Descriptive Statistics --------------------------------------------------------------

# Majority of rulings favor complainant on at least one legal claim
length(d2$prop_adverse[d2$outcome=="Ruling"&d2$prop_adverse>0])/nrow(d2[d2$outcome=="Ruling",])

# Reversed on appeal
length(d2$d_reverse_appeal[d2$d_reverse_appeal==1])/nrow(d2[d2$outcome=="Ruling",])


# Compliance scores
summary(d1$comply_scm)
table(d1$comply_b); 56/120
table(d1$comply_n); 38/120; 58/120

# Third parties
length(d1$thirdparties[d1$thirdparties==0])
length(d1$thirdparties[d1$thirdparties>=1&d1$thirdparties<=5])

# EU dummy
table(d1$is_EU2)

# Federalism and Veto Points correlated
cor(data.frame(d1$federal2, d1$polconiii2, d1$polcon_eu2))

# EU versus Non-EU compliance rates - no significant difference in means
t.test(x=d1$comply_b[d1$is_EU2==1], y = d1$comply_b[d1$is_EU2==0])

# Disputes settled early versus litigated (import restriction disputes)
d.settle = d2[d2$npanel_r==0,]; nrow(d.settle)							# 203 settled cases
d.lit = d2[d2$npanel_r>=1,]; nrow(d.lit)								# 149 litigated cases

t.test(d.settle$polcon_eu2, y = d.lit$polcon_eu2)
t.test(d.settle$polity2, y = d.lit$polity2)



######################################################################################
#### TABLES 1 and 2 --------------------------------------------------------------

## Use PCA to construct domestic constraints score for respondent government (country 2) 

# select variables
pca.vars <-(subset(d1, select = c("ds","federal2","polity2","checks2","polariz2","herfgov2.r","legfralower_eu2")))

# standardize variables
pca.vars$polity2 <- (pca.vars$polity2-mean(pca.vars$polity2, na.rm=T))/sd(pca.vars$polity2, na.rm=T)
pca.vars$federal2 <- (pca.vars$federal2-mean(pca.vars$federal2, na.rm=T))/sd(pca.vars$federal2, na.rm=T)
pca.vars$checks2 <- (pca.vars$checks2-mean(pca.vars$checks2, na.rm=T))/sd(pca.vars$checks2, na.rm=T)
pca.vars$polariz2 <- (pca.vars$polariz2-mean(pca.vars$polariz2, na.rm=T))/sd(pca.vars$polariz2, na.rm=T)
pca.vars$herfgov2.r <- (pca.vars$herfgov2.r-mean(pca.vars$herfgov2.r, na.rm=T))/sd(pca.vars$herfgov2.r, na.rm=T)
pca.vars$legfralower_eu2 <- (pca.vars$legfralower_eu2-mean(pca.vars$legfralower_eu2, na.rm=T))/sd(pca.vars$legfralower_eu2, na.rm=T)

# conduct PCA and attach components 1 and 2 to data frame
pca <- princomp(~polity2+federal2+checks2+polariz2+herfgov2.r+legfralower_eu2,data=pca.vars, na.action=na.exclude)
summary(pca)
d1$pca.c1 <- pca$scores[,2]
d1$pca.c2 <- pca$scores[,1]



## Table 1 -- Probit models	##

r1.1 <-probit(comply_b ~ polconiii2,dat=d1)
r1.2 <-probit(comply_b ~ polconiii2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1, dat=d1)
r1.3 <-probit(comply_b ~ polcon_eu2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1, dat=d1)
r1.4 <-probit(comply_b ~ federal2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1, dat=d1)
r1.5 <-probit(comply_b ~ pca.c1+pca.c2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1, dat=d1)
r1.6 <-probit(comply_b ~ polconiii2+federal2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy, dat=d1[d1$is_EU2==0,])
r1.7 <-probit(comply_b ~ polconiii2+federal2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy, dat=d1[d1$d_flag2==1,])
r1.8 <-probit(comply_b_plty ~ polconiii2+federal2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1, dat=d1)

stargazer(r1.1, r1.2, r1.3, r1.4, r1.5, r1.6, r1.7, r1.8, no.space=T, type="text", digits=2)



## Table 2 -- Multinomial probit models	##

# outcome variable as factor
d1$comply_n <- as.factor(d1$comply_n)
d1$comply_n_plty <- as.factor(d1$comply_n_plty)

r2.1 <-polr(comply_n ~ polconiii2, method="probit", dat=d1)
r2.2 <-polr(comply_n ~ polconiii2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1, method="probit", dat=d1)
r2.3 <-polr(comply_n ~ polcon_eu2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1, method="probit", dat=d1)
r2.4 <-polr(comply_n ~ federal2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1, method="probit", dat=d1)
r2.5 <-polr(comply_n ~ pca.c1+pca.c2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1,method="probit", dat=d1)
r2.6 <-polr(comply_n ~ polconiii2+federal2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy, method="probit", dat=d1[d1$is_EU2==0,])
r2.7 <-polr(comply_n ~ polconiii2+federal2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy, method="probit", dat=d1[d1$d_flag2==1,])
r2.8 <-polr(comply_n_plty ~ polconiii2+federal2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1, method="probit", dat=d1)
r2.9 <- lm(comply_scm ~ polconiii2+federal2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1, data=d1)

stargazer(r2.1, r2.2, r2.3, r2.4, r2.5, r2.6, r2.7, r2.8, r2.9, no.space=T, type="text", digits=2)

cbind(logLik(r2.1),logLik(r2.2),logLik(r2.3),logLik(r2.4),logLik(r2.5),logLik(r2.6),logLik(r2.7),logLik(r2.8))


## Predicted Probability in Table 2.2

# select variables 
d.pred <- subset(d1, select=c("ds","syear","c1_comp","c2_resp","comply_b","comply_n","polconiii2","polcon_eu2",
							"lnthird","prop_adverse","gdp1","gdp2","remedy","leg_b","is_EU2","is_EU_US1"))
r2.2 <-polr(comply_n ~ polconiii2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2, method="probit", dat=d.pred)

# lowest to highest VPs 
vp.range <- quantile(d.pred$polconiii2)[c(1,5)]
temp.data <- data.frame(polconiii2=vp.range, 
						rbind(	apply(X = d.pred[,8:15], MARGIN = 2, FUN = mean),
								apply(X = d.pred[,8:15], MARGIN = 2, FUN = mean)))

							
pred.r2.2 <- predict(r2.2,temp.data,type="p")

diff.2 <- (pred.r2.2[1,3] - pred.r2.2[1,2]) # highest VP, transition 1|2
diff.1 <- (pred.r2.2[2,3] - pred.r2.2[2,2])	# lowest VP, transition 1|2
diff.in.diff <- diff.1-diff.2
diff.in.diff





#### FIGURE 3 --------------------------------------------------------------

# probit model Table 1.3
fig3.fit <-glm(comply_b ~ polcon_eu2+lnthird+prop_adverse+gdp2+gdp1+leg_b+remedy+is_EU2+is_EU_US1,
			family = binomial("probit"), dat=d.pred)

# predicted data
range(d.pred$polcon_eu2)
x.range <- seq(0,1,.005)
temp.data <- data.frame(polcon_eu2=x.range, 
				lnthird=rep(mean(d.pred$lnthird), length(x.range)),
				prop_adverse=rep(mean(d.pred$prop_adverse), length(x.range)),
				gdp2=rep(mean(d.pred$gdp2),length(x.range)),
				gdp1=rep(mean(d.pred$gdp1),length(x.range)),
				remedy=rep(mean(d.pred$remedy),length(x.range)),
				leg_b=rep(mean(d.pred$leg_b),length(x.range)),
				is_EU2=rep(mean(d.pred$is_EU2),length(x.range)),
				is_EU_US1=rep(mean(d.pred$is_EU_US1),length(x.range)))
					
predicted.data <- as.data.frame(predict(fig3.fit, newdata=temp.data, 
					type="link", se=T))
new.data <- cbind(temp.data, predicted.data)

# confidence intervals with jitter

std.99 <- qnorm(0.99 / 2+.5)
std.95 <- qnorm(0.95 / 2+.5)
std.9 <- qnorm(0.9 / 2+.5)

jtr <- rnorm(n=length(x.range),mean=0,sd=.002)
new.data$ymin1 <- fig3.fit$family$linkinv(new.data$fit - std.9 * new.data$se)-jtr
jtr <- rnorm(n=length(x.range),mean=0,sd=.002)
new.data$ymax1 <- fig3.fit$family$linkinv(new.data$fit + std.9 * new.data$se)-jtr
jtr <- rnorm(n=length(x.range),mean=0,sd=.002)
new.data$ymin2 <- fig3.fit$family$linkinv(new.data$fit - std.95 * new.data$se)-jtr
jtr <- rnorm(n=length(x.range),mean=0,sd=.002)
new.data$ymax2 <- fig3.fit$family$linkinv(new.data$fit + std.95 * new.data$se)-jtr
jtr <- rnorm(n=length(x.range),mean=0,sd=.002)
new.data$ymin3 <- fig3.fit$family$linkinv(new.data$fit - std.99 * new.data$se)-jtr
jtr <- rnorm(n=length(x.range),mean=0,sd=.002)
new.data$ymax3 <- fig3.fit$family$linkinv(new.data$fit + std.99 * new.data$se)-jtr

new.data$fit <- fig3.fit$family$linkinv(new.data$fit)  # Rescale to 0-1


# Extract example points  	

ex.pts <- as.data.frame(matrix(NA,6,5))
colnames(ex.pts) <- c("ds","syear","respondent","polcon_eu2","prob_comply")

ex.pts[1,1:4] <- d.pred[23,c(1,2,4,8)] 
ex.pts[1,5] <- new.data[new.data$polcon_eu2==0.210,10]

ex.pts[2,1:4] <- d.pred[92,c(1,2,4,8)] 
ex.pts[2,5] <- new.data[new.data$polcon_eu2==0.290,10]

ex.pts[3,1:4] <- d.pred[120,c(1,2,4,8)] 
ex.pts[3,5] <- new.data[new.data$polcon_eu2==0.400,10]

ex.pts[4,1:4] <- d.pred[66,c(1,2,4,8)] 
ex.pts[4,5] <- new.data[new.data$polcon_eu2==0.635,10]

ex.pts[5,1:4] <- d.pred[98,c(1,2,4,8)] 
ex.pts[5,5] <- new.data[new.data$polcon_eu2==0.685,10]

ex.pts[6,1:4] <- d.pred[119,c(1,2,4,8)] 
ex.pts[6,5] <- new.data[new.data$polcon_eu2==0.835,10]

ex.pts[,4:5] <- round(ex.pts[,4:5], digits=2)


# Plot everything, Figure 3, Predicted Effects

p <- ggplot(d.pred, aes(x= polcon_eu2, y= comply_b)) 
p +  coord_cartesian(xlim = c(.1,.9), ylim = c(.02,.98))+
  	geom_ribbon(data=new.data, aes(y=fit, ymin=ymin1, ymax=ymax1), alpha=0.2, fill="grey30") + 
  	geom_ribbon(data=new.data, aes(y=fit, ymin=ymin2, ymax=ymax2), alpha=0.2, fill="grey30") + 
    geom_ribbon(data=new.data, aes(y=fit, ymin=ymin3, ymax=ymax3), alpha=0.01, fill="grey30") + 
  	geom_line(data=new.data, aes(y=fit), color="grey40",size = 1.1) +
  	labs(x="Veto Points",y="Pr(Comply)",size=5) +

    geom_text(aes(x = ex.pts[1,4], y = ex.pts[1,5], label = "Korea '97"), size=5, nudge_y = 0)+
    geom_text(aes(x = ex.pts[2,4], y = ex.pts[2,5], label = "Mexico '04"),size=5, nudge_y = -.02)+
    geom_text(aes(x = ex.pts[3,4], y = ex.pts[3,5], label = "US '09"), size=5, nudge_y = .02)+
    geom_text(aes(x = ex.pts[4,4], y = ex.pts[4,5], label = "Japan '02"),size=5, nudge_y = .05)+
    geom_text(aes(x = ex.pts[5,4], y = ex.pts[5,5], label = "Brazil '05"),size=5, nudge_x = -.03)+
    geom_text(aes(x = ex.pts[6,4], y = ex.pts[6,5], label = "EU* '09"), size=5, nudge_y = -.02)+

    geom_rug(sides="b",color="black")+
    theme_light()+
    theme(axis.title=element_text(size=14),axis.text=element_text(size=12))








######################################################################################
#### ROBUSTNESS CHECKS --------------------------------------------------------------

## Appendix A4: Nominal Compliance ---------------------------------------------------


## Table A4.1 -- Brewster and Chilton -- Nominal Measure of Compliance in US Cases

d.US <- d2[d2$ccode2==2,]
bc_comply_l <- c(2,4,24,33,58,61,99,108,136,162,138,165,166,177,178,179,192,202,212,213,236,247,257,
				264,277,311,248,249,251,252,253,254,258,259,268,335,343,345,383,392,402)
bc_nocomply_l <- c(160,176,184,217,234,267,285,294,322,344,350,353,379,382,404)
d.US$bc_compliance_1_0 <- ifelse(is.element(d.US$ds, bc_comply_l ), 1, ifelse(is.element(d.US$ds, bc_nocomply_l), 0, NA))

bob <- na.omit(data.frame(d.US$ds, ifelse(d.US$comply_n>0,1,0), d.US$bc_compliance_1_0))
names(bob) <- c("ds", "comply_SCM", "comply_BC")
bob$mismatch <- ifelse(bob$comply_SCM==bob$comply_BC,"","*")
bob[,c(2:3)] <- ifelse(bob[c(2:3)]==1,"comply", "no")
names(bob) <- c("Dispute", "Trade-Based Measure", "Nominal Measure (B&C)","Mismatch?")
print(bob, row.names=F)


## Table A4.2 -- Regression Using B&C Nominal Measure of Compliance, US Only

reg.us1 <- probit(bc_compliance_1_0 ~ leg_b, dat=d.US)
reg.us2 <- probit(bc_compliance_1_0 ~ leg_b+ polconiii2+ prop_adverse, dat=d.US)


# variables and labels
cov_keep <- c(	"polconiii2","prop_adverse","leg_b")				
cov_order <- c("^leg_b$","^polconiii2$","^prop_adverse$")
cov_lab <- c("Legislative","U.S. Veto Points","\\% Adverse")

stargazer(reg.us1, reg.us2, no.space=T, digits=2, type="text",
			keep = cov_keep, order = cov_order, covariate.labels = cov_lab)



## Appendix B1: Selection into Rulings ----------------------------------------------

# Systemic disputes often produce adverse rulings
cor(d2$systemic,d2$adverse)

# Systemic disputes uncorrelated with compliance
cor(d2$systemic,d2$comply_b, use = "complete.obs")


heck1.1 <- heckit(selection = adverse ~ systemic + gdp1 
										+ lnthird + idealpointdistance
										+ elecyr1 + pr1 + W1   
										+ exrate_ratio, 
                 comply_b ~ polconiii2+lnthird+prop_adverse+gdp2+gdp1+remedy+leg_b ,
                 data = d2, method = "2step")

heck2.1 <- heckit(selection = adverse ~ systemic + gdp1  
										+ lnthird + idealpointdistance 
										+ all_case + ag_case 
										+ service_case + sps_case, 
                 comply_b ~ polconiii2 +lnthird+prop_adverse+gdp2+gdp1+remedy+leg_b,
                 data = d2, method = "2step")

# variables and labels
cov_outcome <- c("Constant", "polconiii2", "lnthird", "prop_adverse", "gdp1", "gdp2", "remedy", "leg_b")
cov_o_order <- c("^polconiii2$","^lnthird$","^prop_adverse$","^gdp2$","^gdp1$","^remedy$","^leg_b$","^Constant$")
cov_o_lab <- c("Respondent Veto Points", "Third Parties", "% Adverse Ruling","Respondent GDP", "Complainant GDP", "Remedy", "Legislative Measure", "Constant")

cov_selection <- c("Constant", "systemic", "gdp1", "lnthird", "idealpointdistance", "elecyr1", "pr1", "W1", "exrate_ratio","all_case","ag_case","service_case","sps_case")
cov_s_order <- c("^systemic$","^lnthird$","^idealpointdistance$","^gdp1$","^elecyr1$","^pr1$","^W1$","^exrate_ratio$",
				"^all_case$","^ag_case$","^service_case$","^sps_case$")
cov_s_lab <- c("Systemic Interest", "Third Parties", "Ideal Point Distance", "Complainant GDP", "Complainant Election Year", "Complainant PR", "Complainant MWC", "Exchange Rate Ratio", "All Products","Agriculture","Services", "Sanitary and Phytosanitary")


# Outcome stage table
stargazer(heck1.1, heck2.1, selection.equation=FALSE, no.space=TRUE, type="text",
			dep.var.caption = "Outcome Equation: Did the respondent comply?",
			keep = cov_outcome, order = cov_o_order, covariate.labels = cov_o_lab)   

# Selection stage table
stargazer(heck1.1, heck2.1, selection.equation=TRUE, no.space=TRUE, type="text",
			dep.var.caption = "Selection Equation: Is there an adverse ruling?",
			keep = cov_selection, order = cov_s_order, covariate.labels = cov_s_lab)


## Appendix Table B2: Selection into Litigation ----------------------------------------


# Using data from Davis (2012), compare negotiated to litigated disputes

# Descriptive stats using complete observations for first specification
d.sub <-na.omit(subset(d.davis, select = c("dsu_ts", "polconiii2","lnpolcon","s301_all",
										"impolicy","distortion","progress","duration")))

mean(d.sub$polconiii[d.sub$dsu==1]); length(d.sub$polconiii[d.sub$dsu==1])
mean(d.sub$polconiii[d.sub$dsu==0]); length(d.sub$polconiii[d.sub$dsu==0])
t.test(x = d.sub$polconiii[d.sub$dsu==1], y = d.sub$polconiii[d.sub$dsu==0])


r1 <-probit(dsu_ts ~ polconiii2+
					lnpolcon+s301_all+impolicy+distortion+progress+duration,
					data=d.davis)

r2 <-probit(dsu_ts ~ polconiii2+
					lnpolcon+s301_all+lnprod+lnworldX+mpen_oecd+impolicy+distortion+progress+duration,
					data= d.davis)

r3 <-probit(dsu_ts ~ polconiii2+
					lnpolcon+s301_all+lnprod+lnworldX+mpen_oecd+impolicy+distortion+progress+duration,
					data= d.davis[d.davis$polity2>6,])
summary(r1)
summary(r2)
summary(r3)


# variables and labels
cov_keep <- c("Constant", "polconiii2", "lnpolcon", "lnprod", "s301_all", "lnworldX", "mpen_oecd", "impolicy","distortion","progress","duration")
cov_order <- c(	"^polconiii2$", "^lnpolcon$", "^s301_all$","^lnprod$",
				"^lnworldX$", "^mpen_oecd$", "^impolicy$","^distortion$",
				"^progress$","^duration$","^ Intercept$")
cov_lab <- c(	"Repondent Veto Points", "US Political Contributions", "Section 301","Production Value",
				"Exports Value", "MPEN (partner)", "Import Policy","Distortion",
				"Progress","Duration","Intercept")



stargazer(r1, r2, r3, no.space=T, type="text", dep.var.caption = "File WTO Complaint?",
			keep = cov_keep, order = cov_order, covariate.labels = cov_lab)


## Appendix Table B3: Robustness to Additional Controls ----------------------------------------


reg.b31 <-probit(comply_b ~ polconiii2 
							+lnthird+prop_adverse+gdp2+gdp1+leg_b
							+comp_past_comp+resp_past_resp+nclaims,
							dat=d1)
							
reg.b32 <-probit(comply_b ~ polconiii2 
							+lnthird+prop_adverse+leg_b+is_EU2
							+imports_gdp1+exports_gdp1+exports_gdp2+gdp_ratio,
							dat=d1)
														
reg.b33 <-probit(comply_b ~ polconiii2  
							+lnthird+prop_adverse+gdp2+gdp1+leg_b
							+ elecyr1 + elecyr2,
							dat=d1)

reg.b34 <-polr(as.factor(comply_n) ~ polconiii2
							+lnthird+prop_adverse+gdp2+gdp1+leg_b
							+comp_past_comp+resp_past_resp+nclaims,
							method="probit", 
							dat=d1)


reg.b35 <-polr(as.factor(comply_n) ~ polconiii2
							+lnthird+ prop_adverse+leg_b+is_EU2
							+imports_gdp1+exports_gdp1+exports_gdp2+gdp_ratio,
							method="probit", 
							dat=d1)


reg.b36 <-polr(as.factor(comply_n) ~ polconiii2 
							+lnthird+prop_adverse+gdp2+gdp1+leg_b
							+elecyr1+elecyr2,
							method="probit", 
							dat=d1)


# variables and labels
cov_keep <- c("polconiii2","lnthird","prop_adverse","nclaims"
							,"leg_b","is_EU2"
							,"gdp1", "gdp2", "gdp_ratio" 
							,"comp_past_comp","resp_past_resp" 
							,"imports_gdp1","exports_gdp1","exports_gdp2"  
							,"elecyr1","elecyr2")
cov_order <- c("^polconiii2$","^lnthird$","^prop_adverse$"
							,"^leg_b$","^is_EU2$"
							,"^nclaims$","^comp_past_comp$","^resp_past_resp$" 
							,"^gdp1$", "^gdp2$","^gdp_ratio$" 
							,"^imports_gdp1$","^exports_gdp1$","^exports_gdp2$"  
							,"^elecyr1$","^elecyr2$")
cov_lab <- cov_order <- c("Respondent Veto Points","Third Parties","Adverse Ruling"
							,"Legislative Measure","EU Respondent"
							,"Claims","Comp. Past Disputes","Resp. Past Disputes" 
							,"Complainant GDP", "Respondent GDP","GDP Ratio" 
							,"Complainant Imports","Complainant Exports","Respondent Exports"  
							,"Comp. Election Year","Resp. Election Year")




stargazer(reg.b31, reg.b32, reg.b33, reg.b34, reg.b35, reg.b36, no.space=T, type = "text", 
			keep = cov_keep, order = cov_order, covariate.labels = cov_lab)


logLik(reg.b34)
logLik(reg.b35)
logLik(reg.b36)



## Appendix Table B4: Groups of WTO Disputes ----------------------------------------

## Sandwich estimator cluster robust standard errors

clx.glm <- function(dat, l.model, group.id, df.adj=0){
				bread <- vcov(l.model)
				est.fun <-estfun(l.model)
				meat <- t(est.fun)%*%est.fun
				sandwich <- bread%*%meat%*%bread
				coeftest(l.model, sandwich)
				fc <- group.id
				m <- length(unique(fc))
				k <- length(coef(l.model)) + df.adj
				u <- estfun(l.model)
				u.clust <- matrix(NA, nrow=m, ncol=k)
				for(j in 1:k){
					u.clust[,j] <- tapply(u[,j],fc,sum)
					}
				cl.vcov <- vcov(l.model) %*% ((m/(m-1))*t(u.clust)%*% u.clust) %*% 	vcov(l.model)
				r.cl <- coeftest(l.model, cl.vcov)
				return(r.cl)		
			} # end clx.glm()


d.sub <-(subset(d1, select = c("ds","idx", "comply_b", "comply_n", "polconiii2", "lnthird", "prop_adverse", "gdp2","gdp1","remedy", "leg_b","is_EU2")))

reg.b41 <-probit(comply_b ~ polconiii2 
							+lnthird+prop_adverse+gdp1+gdp2+leg_b+remedy+is_EU2,
							dat=d.sub)
clx.b41 <- clx.glm(dat=d.sub, l.model=reg.b41, group.id=d.sub$idx)

reg.b42 <-polr(as.factor(comply_n) ~ polconiii2
							+lnthird+prop_adverse+gdp1+gdp2+leg_b+remedy+is_EU2,
							method="probit", 
							dat= d.sub)							
clx.b42 <- clx.glm(dat=d.sub, l.model=reg.b42, group.id=d.sub$idx, df.adj=2)


# variables and labels
cov_keep <- c("polconiii2","lnthird","prop_adverse","nclaims"
							,"leg_b","is_EU2","remedy"
							,"gdp1", "gdp2")
cov_order <- c("^polconiii2$","^lnthird$","^prop_adverse$"
							,"^gdp1$", "^gdp2$"
							,"^remedy$","^leg_b$","^is_EU2$")
cov_lab  <- c("Respondent Veto Points","Third Parties","Adverse Ruling"
							,"Complainant GDP", "Respondent GDP"
							,"Remedy","Legislative Measure","EU Respondent")
										
stargazer(clx.b41, clx.b42, no.space=T,  
			keep = cov_keep, order = cov_order, covariate.labels = cov_lab)



## Appendix Table B5: Single-step Linear Regression ----------------------------------------

# full bilateral product-level trade data set by dispute-partner-year
nrow(d3)
length(unique(d3$ds))
length(unique(d3$year))


## Clustered SE function - Linear
clx   <- function(dat,fm, cluster){
  attach(dat, warn.conflicts = F)
  library(sandwich)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
  

## Subset data by control variables used

keep_vars <- c("ds","year","idx","ccode_resp","share","polconiii2","polconiii_resp","resp_postT","postT","resp","gdp1_n","gdp2_n",
				"agriculture_va2","industry_va2","manufacturing_va2","unemployment2","polity2")
						
d3.1 <-  na.omit(subset(d3,select = keep_vars))		
				
d3.2 <-  na.omit(subset(d3,select = c(keep_vars,
				"gdp_growth2","gdppc2","trade2")))		

d3.3 <-  na.omit(subset(d3,select = c(keep_vars,
				"lnthird","prop_adverse","leg_b","is_EU2")))

d3.4 <-  na.omit(subset(d3,select = c(keep_vars,
				"lnthird","prop_adverse","leg_b","is_EU2","is_EU_US1",
				"AGcase","SPScase")))


# Reduced form with basic SCM controls				
r1 <-lm(share ~ polconiii_resp*resp_postT +
				gdp1_n+gdp2_n+
				agriculture_va2+industry_va2+manufacturing_va2+unemployment2,
				dat=d3.1)
r1.c <- clx(dat=d3.1, fm = r1, cluster = ds)		
		
# Reduced form with all the SCM controls
r2 <-lm(share ~ polconiii_resp*resp_postT+
				gdp1_n+gdp2_n+
				agriculture_va2+industry_va2+manufacturing_va2+unemployment2+
				gdp_growth2 + I(log(gdppc2)) + I(log(trade2)) + polity2,
				dat=d3.2)
r2.c <- clx(dat=d3.2, fm = r2, cluster = ds)		

# Full regression with essential dispute controls and  SCM controls
r3 <-lm(share ~ polconiii_resp*resp_postT+
				lnthird*resp_postT+ prop_adverse*resp_postT+leg_b*resp_postT+is_EU2*resp_postT+
				gdp1_n+gdp2_n+
				agriculture_va2+industry_va2+manufacturing_va2+unemployment2+polity2,
				dat=d3.3)
r3.c <- clx(dat=d3.3, fm = r3, cluster = ds)		

# Full regression with extra dispute controls and  SCM controls
r4 <-lm(share ~ polconiii_resp*resp_postT+
				lnthird*resp_postT+ prop_adverse*resp_postT+leg_b*resp_postT+is_EU2*resp_postT+
				AGcase*resp_postT+SPScase*resp_postT+
				gdp1_n+gdp2_n+
				agriculture_va2+industry_va2+manufacturing_va2+unemployment2+polity2,
				dat=d3.4)
r4.c <- clx(dat=d3.4, fm = r4, cluster = ds)		



# variables and labels
cov_keep <- c(	"polconiii_resp","resp_postT", "gdp1_n","gdp2_n","gdp_growth2",
				"gdppc2","trade2","agriculture_va2","industry_va2","manufacturing_va2",
				"unemployment2","polity2","lnthird","prop_adverse","leg_b",
				"is_EU2","is_EU_US2","remedy","AGcase","SPScase")
				
cov_order <- c("^polconiii_resp:resp_postT$","^polconiii_resp$","^resp_postT$", 
				"^resp_postT:lnthird$","^lnthird$",
				"^resp_postT:prop_adverse$","^prop_adverse $",
				"^resp_postT:leg_b$","^leg_b$",
				"^resp_postT: is_EU2 $","^ is_EU2 $",
				"^resp_postT:AGcase$","^AGcase$",
				"^resp_postT:SPScase$","^SPScase$", 
				"^gdp1_n$","^gdp2_n$",
				"^agriculture_va2$","^industry_va2$","^manufacturing_va2$",
				"^unemployment2$","^polity2$",
				"^gdp_growth2$","^I(log(gdppc2))$","^I(log(trade2))$")

cov_lab <- c(	"Respondent Veto Points$\\times$Treated","Respondent Veto Points", "Treated",
							"Third Parties$\\times$Treated","Third Parties",
							"Adverse Ruling$\\times$Treated","Adverse Ruling",
							"Legislative Measure$\\times$Treated","Legislative Measure",
							"EU Respondent$\\times$Treated","EU Respondent",
							"Agriculture$\\times$Treated","Agriculture",
							"Sanitary $\\&$ Phytosan$\\times$Treated", "Sanitary $\\&$ Phytosan",
							"Complainant GDP", "Respondent GDP",
							"Partner agriculture","Partner industry","Partner manufacturing",
							"Partner unemployment","Partner Polity score", "Partner GDP growth",
							"Partner GDP per capita","Partner trade dependence")



stargazer(r1.c, r2.c, r3.c, r4.c, no.space=T,  
			keep = cov_keep, order = cov_order, covariate.labels = cov_lab)
  

nrow(d3.1); length(unique(d3.1$ds)); length(unique(d3.1$year)); summary(r1)[8] 
nrow(d3.2); length(unique(d3.2$ds)); length(unique(d3.2$year)); summary(r2)[8] 
nrow(d3.3); length(unique(d3.3$ds)); length(unique(d3.3$year)); summary(r3)[8] 
nrow(d3.4); length(unique(d3.4$ds)); length(unique(d3.4$year)); summary(r4)[8] 






## Appendix Table B6.1: Principal Component Analysis ----------------------------------------

# Variable Loadings and Importance of Components
# Relies on main analysis code above
loadings(pca)



## Appendix Figure B7: Veto Players in European Union ----------------------------------------

## EU robustness plot

d.i <- d.sub.eu <- subset(d1, select = c("comply_b","polconiii2","lnthird","prop_adverse","gdp1","gdp2","leg_b","remedy","is_EU2"))


VP <- seq(.05,.95,.05)
stan <- cbind(VP,matrix(NA,length(VP),3))
colnames(stan) <- c("VP","beta","ci5", "ci95")

for(i in 1:length(VP)){
	d.i$VP_i <- ifelse(d.i$is_EU2==1, VP[i], d.i$polconiii2)
	reg_i <- glm(comply_b ~ VP_i 
			+lnthird+prop_adverse+gdp1+gdp2+leg_b+remedy,
			dat=d.i, family = binomial(link="probit"))

	stan[i,2] <- coef(reg_i)[2]
	stan[i,3:4] <- confint(reg_i, level=.9)[2,]
}

stan<-as.data.frame(stan)

## Plot Figure B7 - EU veto points robustness
par(mar=c(6,6,3,3))
plot(stan$VP, stan$beta, "l", ylim = c(-5,2), lwd=2, xlim = c(.05, .95),
			xlab = "Hypothetical Veto Points Score for European Union",
			ylab = "Beta Coeffient on Veto Points, Probit Model with All Controls")
abline(h = 0, lwd=1, col="darkgrey")
lines(x=stan$VP, y=stan$ci5)
lines(x=stan$VP, y=stan$ci95)
abline(v=.2, col = "black", lty=2)



## END











